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Abstract. In [H], we introduced and analyzed an improved Zienkiewicz-Zhu (ZZ) estima¬ 
tor for the conforming linear finite element approximation to elliptic interface problems. 
The estimator is based on the piecewise “constant” flux recovery in the ir(div;0) con¬ 
forming finite element space. This paper extends the results of [5] to diffusion problems 
with full diffusion tensor and to the flux recovery both in piecewise constant and piecewise 
linear H{div) space. 


1 Introduction 

A posteriori error estimation for finite element methods has been extensively studied for the 
past four decades (see, e.g., books by Ainsworth and Oden [I], Babuska and Strouboulis [3], 
and Verfiirth [2S] and references therein). Due to easy implementation, generality, and ability 
to produce quite accurate estimation, the Zienkiewicz-Zhu (ZZ) recovery-based error estimator 
m has been widely adapted in engineering practice and has been the subject of mathematical 
study (e.g., |llEllHl[IIl[I31[IHl[lil231l211l2Hl[2Sl[2Sll3niE2). By first recovering a gradient 

(flux) in the conforming linear vector finite element space from the numerical gradient 
(flux), the ZZ estimator is defined as the norm of the difference between the recovered and 
the numerical gradients/fluxes. 

Despite popularity of the ZZ estimator, it is also well known (see, e.g., IZQl EH) that 
adaptive mesh refinement (AMR) algorithms using the ZZ estimator are not efficient to reduce 
global error for non-smooth problems, e.g., interface problems. This is because they over- 
refine regions where there are only small errors. By exploring the mathematical structure of 
the underlying problem and the characteristics of finite element approximations, in mm we 
identified that this failure of the ZZ estimator is caused by using a continuous function (the 
recovered gradient (flux)) to approximate a discontinuous one (the true gradient (flux)) in the 
recovery procedure. Therefore, to fix this structural failure, we should recover the gradient 
(flux) in proper finite element spaces. More specifically, for the conforming linear finite element 
approximation to the interface problem we recovered the flux in the ih(div; D) conforming finite 
element space. It was shown in [5] that the resulting implicit and explicit error estimators are 
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not only reliable but also efficient. Moreover, the estimators are robust with respect to the 
jump of the coefficients. 

In [8], the implicit error estimator requires solution of a global minimization problem, 
and the explicit error estimator uses a simple edge average. This averaging approach of the 
explicit estimator is limited to the Raviart-Thomas {RT) [7] element of the lowest order, i.e., 
the piecewise “constant” vector. In this paper, we introduce a general approach of constructing 
explicit flux recoveries using either the piecewise “constant” vector or the piecewise linear vector 
(the Brezzi-Douglas-Marini (BDM) [7] element of the lowest order) for the diffusion problem 
with the full diffusion coefficient tensor. With the recovered fluxes, the improved ZZ estimators 
are defined as a weighted norm of the difference between the recovered and the numerical 
fluxes. These estimators are theoretically shown to be locally efficient and globally reliable. 
Moreover, when the diffusion coefficient is piecewise constant scalar and its distribution is 
locally quasi-monotone, these estimators are robust with respect to the size of jumps. For 
a benchmark test problem, whose coefficient is not locally quasi-monotone, numerical results 
also show the robustness of the estimators. 

The paper is organized as follows. Section 2 describes the diffusion problem, variational 
form, and conforming finite element approximation. Section 3 describes the a posteriori error 
estimators of the ZZ type and two counter-examples. Two explicit flux recoveries and their 
corresponding improved ZZ estimators are introduced in section 4. Efficiency and reliability 
of those estimators are established in section 5. Section 6 is devoted to explicit formulas of 
the recovered fluxes, the indicators, and the estimators. Finally, section 7 presents numerical 
results on the Kellogg’s benchmark test problem. 

2 Finite Element Approximation to Diffusion Problem 

Let H be a bounded polygonal domain in 51?'^ with d = 2 or 3, with boundary 90 = F^ U F^, 
F^ n Fj^ = 0, and measrf_i (F^) ^ 0, and let n be the outward unit vector normal to the 
boundary. Consider diffusion equation 

-V • {A{x)Vu) = f in O (2.1) 


with boundary conditions 


—AVu • n = oil and u = g^ on F^, (2.2) 

where the V- and V are the divergence and gradient operators, respectively, and / G L^(II). 
In this paper, we consider only simplicial elements. Let T = {K} be a regular triangulation of 
the domain Q, and denote by the diameter of the element K. For simplicity of presentation, 
assume that and g^^ are piecewise affine functions and constants, respectively, and that A 
is a piecewise constant matrix that is symmetric, positive definite. 

Let 

= {v : v = g^ on F^} and i7^(II) = 

Then the corresponding variational problem is to find u G L7j^(II) such that 

a{u, v) = {AVu,Vv) = {f,v) - {g^,v)r^ = f{v) VuGiL^(O), (2.3) 

where (•, •)tj is the inner product on the domain uj. The subscript oj is omitted when a; = H. 
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For each K £ T, let Pk{K) be the space of polynomials of degree less than or equal to k. 
Denote the linear conforming finite element space m associated with the triangulation T by 


S = {v £ H^{n) : v\k e Pi{K) yK£r}. 


Let 


■’ 9,0 


= {v £ S : V = on T and 5^ = S^ 


0,0 


Then the conforming finite element approximation is to seek Uj- £ such that 

V^;) = f{v) 


yv £ s^. 


(2.4) 


3 ZZ Error Estimator and Counter Examples 


. Denote 

the numerical gradient and the numerical flux by 


Let Uj. £ be an approximation to the finite element solution Uj- £ Sg^^ of (2.4 


and = —AVu^, (3-1) 

respectively, which are piecewise constant vectors in with respect to the triangulation T. 
It is common in engineering practice to smooth the piecewise constant gradient or flux in a 
post-process so that they are continuous. More precisely, denote by M the set of all vertices 
of the triangulation T. For each vertex z £ M, denote by Uz the union of elements having z 
as the common vertex. The recovered (smoothed) gradient or flux are defined as follows: for 
r = or a-,j- 

G{t) £ S'^ C C'°(D)'^ with nodal values G{t){z) = -r^ f t dx M z £ M, (3.2) 

I \ j OJz 

where |a;^| = measd(w^). There are many post-processing, recovery techniques (see survey 
article |29] by Zhang and references therein). With the recovered gradient (flux), the ZZ error 
indicator and estimator are defined as follows: 


^zz,K = l|G'('^) - '^\\o,K y K £T and = ||G(t) - t||o,q 
for T = oi respectively. 

Despite many attractive features of the ZZ error estimator, it is well known (see Figures 
in section 7 that adaptive mesh refinement (AMR) algorithms using the above ZZ estimator 
are not efficient to reduce global error for interface problems. In this section, we demonstrate 
this failure of the ZZ estimator by two counter-examples for which the finite element solution 
is exact while the ZZ indicators along interface could be arbitrarily large. 

The first example is a one-dimensional interface problem defined on the domain D = (0, 1) 
with the Dirichlet boundary condition: 

tt(0) = 0 and rt(l) = (A:-|-l)/2 

for an arbitrary constant A: > 1 and with piecewise constant diffusion coefficient 
A = 1 in (0, 1/2) and A = A: > 0 in (1/2, 1). 
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The exact solution and its derivative of this example are piecewise linear and piecewise constant 
functions, respectively, depicted in Figures 1 and 2 and given by 

kx X G (0, 1/2], 

u = < k — 1 ~ 


X + 


XG (1/2, 1), 


k 

1 


X G (0, 1/2), 
xG(1/2, 1). 




Figure 1: solution u = 


Figure 2: u'^{x) and G{u'^) 


For any triangulation T with x = 1/2 as one of its vertices, the conforming linear finite 
element approximation is identical to the exact solution: = u, and hence the true error 

equals to zero. Without loss of generality, assume that the size of two interface elements is h. 
Then the recovered gradient is depicted in Figure 2 and its value at x = 1/2 is {k + l)/2. A 
simple calculation yields the ZZ error indicator: 


= \\G{u' 


1 


T' 


- < 


2^/3 

1 

2\/3 

0 , 


{k-l)h^G^ 

{k-l)h^G^ 


xG (1/2-/i, 1/2), 

X G (1/2, 1/2+ /i), 
otherwise. 


Hence, no matter how small the mesh size h is, the ZZ indicators at two interface elements 
could be arbitrarily large. 

For this simple one-dimensional example, to overcome the inefficiency of the estimator, 
one may use the ZZ estimator based on the fiux. However, this idea could not be extended 
to two or three dimensions. To see this, consider the second example defined on the domain 
H = (—1, 1)^ with scalar piecewise constant diffusion coefficient 


A = kl for y > 0 and A = I for y < 0, 


where fc > 1 is an arbitrary constant. Choose proper Dirichlet boundary data such that the 
exact solution of (2.11 is piecewise linear function given by 


X -I- y if y > 0, 

X + ky if y < 0. 


The conforming linear finite element approximation on any triangulation aligned with the 
interface y = 0 is identical to the exact solution, and hence the true error vanishes. Since the 
exact gradient and flux 

r (1, 1)*, for y > 0, f (^’ 

Vu = < and cr = —AVu = < 

( (1, A:)*, for y < 0 ( (1, A:)*, for y < 0 

are not continuously across the interface, similar calculation to the first example yields that 
the ZZ error indicators on the interface elements based on continuous gradient or flux recovery 
could be arbitrarily large no matter how small the mesh sizes of the interface elements are. 
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4 Improved ZZ Estimators 


The second example of the previous section shows that the gradient and the flux of the ex¬ 
act solution of (2.1) are not continuously across interfaces. This means that inefficiency of 
the ZZ estimator is caused by using continuous functions (the recovered gradient or fiux) to 
approximate discontinuous functions (the true gradient or flux). In this section, we introduce 
improved ZZ estimator that is efficient. 

To this end, let 


H{div, n) = {T€ : V • r E t 

with the norm ||T||j^(div;n) = -h ||V • ^ and let 

Q) = {r £ H{div, Q) : r • n|r^ = g^}. 

Denote the ^f(div; D) conforming Raviart-Thomas {RT) and Brezzi-Douglas-Marini {BDM) 
spaces [7] of the lowest order by 

RT = {r E R(div; D) : t\k G RT{K) VRT E T} 

and BDM = {r E //(div; D) : t\k G BDM{K) V RT G T}, 

respectively, where RT{K) = Pq{KY + (a^i, XdY Pq{K) and BDM{K) = Pi{KY. Let 


RTg,N = RT n iLg^ 7 v(div; D) and BDMg^N = BDM n Hg^N{div, D). 


Let u G R^(D) be the exact solution of (2.11, it is well known that the tangential compo¬ 
nents of the gradient and the normal component of the flux are continuous. Mathematically, 
we have 

“u G Vu G R(curl; D)” and <t = —AVu G LI(div; D), 


where R(curl;D) ^ is the collection of vector-valued functions that are square in- 

tegrable and whose curl are also square integrable. This suggests that proper finite element 
spaces for recovering the gradient and the flux are the respective iL(curl;D) and R(div;D) 
conforming finite element spaces. 

For the conforming finite element approximation, the numerical gradient is already in 
iL(curl; D) and, hence, the resulting improved ZZ estimator based on the gradient recovery 
is identical to zero. Since the numerical flux is not in iLg^ 7 v(div; D), the improved ZZ 
estimators introduced in [5] are based on either explicit or implicit flux recoveries in RTg^i^ 
and BDMg^]\f. The explicit recovery is limited to the scalar diffusion coefficient and the RT 
element. The implicit recovery requires to solve the following global minimization problem: 
find E V such that 


l|4 ^/^(crr-d-r)||oR = min ||4 ^/^(r - o-^)|| o,q. 


(4.1) 


where V = RTg^jsr or BDMg^jsf. With the recovered fiux cr^ E V, the improved ZZ estimator 
introduced in [8] is given by 

£(.^r) = M-''Vr-<^T)lkn- (4-2) 
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Even though the solution of (4.11 may be computed efficiently by a simple iterative solver, in 


the remainder of this section, we derive two new explicit and efficient fiux recoveries applicable 
to the problem with the full diffusion tensor based on the respective RT and BDM elements 
of the lowest order. 

Here, we introduce some notations. Denote the set of all edges (faces) of the triangulation 
by <5 := SjUS^USj^, where is the set of all interior element edges (faces), and <5^ and are 
the sets of all boundary edges (faces) belonging to the respective and For each F G £, 

denote by the unit vector normal to F. For each F G £j, let K~ and be two elements 

sharing the common edge (face) F such that the unit outward normal vector of K~ coincides 
with n^; and let and be the vertex of Kp and Kp opposite to F, respectively. When 
F G £j,U £^, is the unit outward vector normal to d£l and denote by the element 
having the edge (face) F and by x](^ the vertex in Kp opposite to F. 

Let 5 be the Kronecker delta function. For each F G £, denote by (f)^ the global nodal 
basis function of RT associated with F, i.e., 


jp, ^ F' g£] 

denote by -0^^, ..., 0^^ global basis functions of BDM satisfying 

0 +---+0 = d) , M F G £\ 

and let 

RTp = span{0^} and FDM^ = span{0^ ..., 0^^}. 
Since is piecewise constant, for any r G RTq^n or BDMg^N, we have that 

-= E -f+E^., ^|F|0^ with Tp G RTp or BDMp, 


(4.3) 


(4.4) 


(4.5) 


Sj U S 


Fes, 


where g^^p = g^^lp and |F| = measd_i (F). 

For any K G F, restriction of the numerical flux a-p on LC is a constant vector and has the 
following representation in RT{K) (see Lemma 4.4 of [5]): 

^t\k ~ E! ^F,k\R\^ pi 
Fc.dK 

where cTp^ = [ap\K ' ^p)p is the normal component of on F. On each interior edge 

(face) F G £j, the normal component of the numerical flux has two values 


= a and at = a 

P PKp P P’Pp 


Denote by 0^ and cj)'^ the restriction of 0^ on and K^, respectively. Then the numerical 
flux also has the following edge (face) representation: 


v = E 

FcE 


a-p with ap = 


d-|F|0-+d+|F|0+, VFeT„ 




vfet^ut^. 


(4.6) 


For any r E RTg^N or BDMg^N, (|4.5[) and (4.6) give 


= E -^f) + E {'^F-^p\R\(t>p) + E {9n,f-^ p) \R\^p^ 

FcEj “ ■ “ ■ 


FCEr 


FCE, 
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which, together with the triangle inequality and the choice of Tp = for all F € 

implies 


^(crr)=mm p (r - ||o,n 

^ H - ^f) llo,a;^ + I] {9n,f “ ^f) \F\^f\\o,K-^ (4-7) 

FgSj^f^^f 


where iVp is the union of elements sharing the edge (face) F for all F £ V = RTg^jsr or 
BDMg^Ff and = RTp or BDMp. 

For each F ^ £j, let a-p G Vf = RTp or BDMp be the solution of the following local 
minimization problem: 


1^ {a-p -ap)\\o^^^ = P ^/^(T-d-p) 


0,CtJr 


(4.8) 


by (4.7), it is then natural to introduce the following edge (face) based estimator and indicators 

i^F - ^f)\\o,u1p FgSi, 

\F\(t>p\\,^^- F€£n, 


I = ^F '«"ith ^p = 

F^£i U £n 


which satisfies 


per ) < 


To introduce the element based estimator, define the recovered flux G RTg^isf or BDMg^F! 
as follows: 

= Y.^f+ H ^p\F\4)p+ Y. 9n,f\T\4>p- 

F&Ej 


(4.9) 

(4.10) 
JMg^N 

(4.11) 


Then the element based indicators and estimator are given by 

- ^r)llo,i^, yKeT and ^ - 


(4.12) 


The minimization problem in (4.8) is equivalent to the following variational problem: find 


a-p G such that 


(a ^&p, r) = (a ^ap, r) 

\ / \ / L 


V r G 


(4.13) 


The local problem in (4.13) has only one unknown if = RTp and d unknowns if = BDMp. 
The explicit formula of the solution a-p will be given in section 6. 


5 Efficiency and Reliability 


This section establishes efficiency and reliability bounds of the indicators and estimators defined 


in (4.9) and (4.12), respectively, for the diffusion problem with the coefficient matrix A being 


locally similar to the identity matrix. 

To this end, for each K £ T, denote by 
eigenvalues of Ak = A\^, respectively. Let 


and \rmn,K the maximal and minimal 


Amax — max Amax.i^ 
K&T 


and 


Amin = min Amin.fC- 
K&T 
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Assume that each local matrix Ak is similar to the identity matrix in the sense that its maximal 
and minimal eigenvalues are almost of the same size, i.e., there exists a moderate size constant 
K > 0 such that 

Amax,i^ < K, V a: e r. ( 5 . 1 ) 


A 


min.K 


Nevertheless, the ratio of the global maximal and minimal eigenvalues, Amax/Amiiu could be 
very large. In order to show that the reliability constant is independent of the ratio, we assume 
that the distribution of AminC^;) is quasi-monotone (see [22] )• 

Let be the set of all interfaces of the diffusion coefficient that are assumed to be aligned 

1 f 

with element interfaces, and denote by fz =-—r / f dx the average of / over Uz- Let 

measrf(wz Jui, 


Hf = 


E 


measrf(a;2) 


1/2 


A 


II/-/. 


min,aj2 


^\\0,uJz + E E \ ^ \\f\\o,K 

zeMniT mm.ii: 


Note that A is a constant matrix in if z £ AA \ (L^ U L/)). 

Remark 5.1. For various lower order finite element approximations, the first term in Hf is 
of higher order than rjp {defined below in (5.5)) for f G L?‘{Fl) and so is the seeond term for 
f G L'P{Q) with p > 2 {see (Hj). 

Theorem 5.2. (Global Reliability) Assume that the distribution of Xmin,K is quasi-monotone. 
Then the error estimators ^ and defined in (4.9) and (4.12) , respeetively, satisfies the fol¬ 
lowing global reliability bound: 


\\A^/^V{u-u^)\\o,n<C + 


and 


\\A^/^V{u-u^)\\o,n<C (i + Hf), 
where the eonstant C depends on the shape regularity of T and k, but not on Amax/Ar 


(5.2) 

(5.3) 


Proof. It follows from (4.6), (4.11), and Young’s inequality that 

2 

{2= ^ A-'/V,-. 


FgSjUSn 


o,n 


( 


= E E (aI -d-^),(d-^, -d-^,))^_+E (aI 

F^EiWEiAFt^QK- 


F<^Si 


F'CdKl 


’ F)\\o,i.up — ( 2 + ^ ^ 



Now, it suffices to prove the validity of ( |5.2[ ). Note that for any K and for any vector filed 
r, we have that 

ikiio,i^ < iiAi^/^T|io,i^ < AEx,i^ iit||o,k, 
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and that 


>-1/2 II II , II .-1/2 II ^ \-l/2 II II 

-^max.iC lpllo,il' < 11^ '^\\0,K < ^ra\n,K lpllo,ii'- 


(5.4) 


With the above inequalities, (5.2) may be proved in a similar fashion as in mm with the 
constant C also depending on k. □ 


In the remaining part of this section, we will establish the efficiency of the indicators 
and given in (4.9) and (4.12), respectively, by proving that the indicators and are 
bounded above by the classical residual based indicators of the flux jump on edges (faces) given 
in (5.5), which is well known to be efficient for interface problems, (i.e., A = a{x) I with a{x) 
being a piecewise constant function). More specifically, Petzoldt (see (5.7) in [22]) proved that 
the edge (face) flux indicator 


Vf = 


1^1 l^f -^JI/'\/«F +"f’ 

Fg8„ 


W 

1 

1 

L' 

F G Sn, 

(5.5) 

0, 

F G Sd, 



where = oi^±^ is locally efficient without assumption on the distribution of the coefficient 
a. More specifically, there exists a constant C > 0 independent of a and the mesh size such 
that 




< C \\a-^/‘^V{u- 


hlr 


'■rhlO.a; 


+ E 4^11/-/- 


K C u)f ^ 


rllo,iC 


(5.6) 


where fj- is the projection of / onto the space of piecewise constant with respect to T. 


Remark 5.3. For the diffusion problem, define t he ed ge (faee) estimator rjp aeeording to (5.5) 
with then the loeal efifieieney in (5.6) holds with a = A and the eonstant C 

depending on k. 


Theorem 5.4. (Local Efficiency) The loeal edge (faee) and element indicators defined in (4.9) 
and (4.12), respectively, are efficient, i.e., there exists a constant C > 0 depending only on the 
shape regularity of T and k such that 


tbdm 


<eJ<c\\\A^/^VeJo,.,+ 


E -^ 11 /-/- 


1/2N 


(y 

\K'CLOp K' 


t\\1,K' 


yF gs 


(5.7) 


and that 


h^ 


1/2N 




°‘k' 




yKGT, (5.8) 


where is the union of all elements that shares at least one edge (face) with K. 


Proof. It follows from (4.6), (4.11), and the triangle inequality that 

= ll^■^^^(^r-^r)lloW < E -^F)llo,i^ < E ViLer. 

FcdK FcdK 
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Hence, (5.81 is a direct consequence of (5.7). 


To prove the validity of (5.7), first note that the first inequality is a direct consequence 
of the minimization problem in (4.8) and the fact that C To prove the second 


inequality in (5.7), without loss of generality, assume that F £ £j and that ^+. 


By (5.4) and the fact that < C with constant C > 0 depending only on the shape 

regularity of T, we have 


C 11^ < P 


tgrt, 


''fIIIO,..;, 


{i; - i;) A-^/^-\F\,l> 


fuO,K- < 


a — a 

F F 


^-,,2 


icin'/-, 


min,iCj^ I"" III 


< C\F\ 


a — a 

F F 




- 1/2 


Combining with Remark 5.3 implies the second inequality in (5.7) and, hence, (5.8). This 
completes the proof of the theorem. □ 


6 Explicit Formulas 


This section presents explicit formulas of the recovered fluxes defined in (4.11) (see (6.1) and 


(6.5)) and the corresponding indicators and estimators defined in (4.9) and ( 4.12[ ), respectively. 
In particular, the explicit formulas for the indicators and, hence, the estimators are written in 
terms of the current approximation Uj- and geometrical information of elements. For simplicity, 
we only consider the two-dimensional case. 

For each edge F € £, denote by and the globally fixed initial and terminal points of 
F, respectively, such that = |F| with being a unit vector tangent 

to F; by p, —t^ p) a unit vector normal to F; and by the opposite vertices of F in 

respectively. 

Denote by Ag^ and Ae^ the nodal basis functions of the continuous linear element associated 
with vertices and of A7, respectively. For any v G denote the formal adjoint of 

the curl operator by V"*"?; = {dv/dy,—dv/dxY. For the RT space of the lowest index, the 
nodal basis function associated with F G £ is given by 




For the BDM space of the lowest index, two basis functions associated with the edge F G £ 
are given by 

-^s,p = ^Sp^^Kp and -0^^^ =-Ae^V^Ag^, 
respectively, which satisfy 




= Ag^5^^,/|F'| and (0^^^ • n^/) 


= AeP 


p FF' I 


IF'I 


for any F' G £. It is easy to check that (4.3) and (4.4) hold 


A detaied MATLAB implementation of BDM and RT mixed finite element methods can 
be found in m- 
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6.1 Indicator and Estimator Based on RT 


For all F £ let 


7 ^ = {a and 


7, 


+7J 

Using the basis function cj) defined above, a straightforward calculation gives 


7f = 


48| iiT; 


/ 

i: ||a-‘/^(x±-x±) 

\F'CdK^ 


+ 


71 - 1/2 


E x^,-3x^ 

\F'CdK^ 


where || • || is the standard Euclidean norm in TZ^. Solving the local problem in (4.13) with 
Vp = RTp gives the following recovered flux in RTg^^: 


»r= E E w+ E s„.,\F\4,; 

F&Sj 


( 6 . 1 ) 


with the normal component of the recovered flux, d))*, on each edge F £ £j given by the 
following weighted average: 

The edge indicator has the following explicit formula: 

7/2 


( 6 . 2 ) 


C = { 0 


-^^1 1^1 ((1 -aF)^7 F+«f7J) > F€£i, 

F e £d, 

-5jv,fI i^i yE’ ^ ^ 


Next, we introduce explicit formula of the element indicator in terms of the current 
approximation Up and geometrical information of elements. To this end, for any K € F, 
denote the sign function by 


sign (F) = 


1 if n^=n^|^. 


-1 if n^ = -n^|^. 


VF c dK, 


where is the unit outward vector normal to K, and let 


^FF' = E ^F" - E ^F» - + E (^F- - - X^), 

\F’'cdK j \F"cdK j F"cdK 

Tpp, = nM E Xf"-3x^)> and 

\F"CdK j 

For each element K the indicator is given by 


C' = (7“''»/. + 2 + {AVu^, v«o,) 


1/2 


(6.3) 
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with explicit formula for each term as follows 


(a-' 


48\K\ 


sign;^(F)sign^(F') |-F||F'|<t;* 


FcdK F'cdK 


(<, 


12\K\ 


1 


J2 sign^(F)sign^(F') 


FF' ’ 


FcdK F'cdK 


and {AVu^,Vu^)^ =-— J2 J2 ^r(xF')signi^(F)sign;^(F') |F||F'|S'^^,. 

^ I FcdK F'cdK 


Remark 6.1. For interface problems, the recovered flux in (6.11 and the resulting estimator 


defined in (4.121 are equivalent to those introduced and analyzed in [S]. To see that, let A\^ = 
aj^I for any K G T, where aj^ and I are constant and the identity matrix, respectively. Let 


then 


a = a _ and a~f = a , , 


1 1 


For a regular triangulation, the ratio of {flpp, o,nd (fpp, 4>p^ Is bounded above and 

below by constants. Thus 


% 


Oil 


and 1 — a„ = 


V 


ct. 


(6.4) 


^F + “f + “f ^F +^F "f + 

[Here, we use x ^ y to mean that there exist two positive constants Ci and C 2 independent 


of the mesh size such that Cix < y < C 2 X.) (6.4| indicates that the weights in (6.1) may 

a~f a~ 

be replaced by the respective _ —— and ——-—— Hence, it is equivalent to the explicit 


ap+a-f 




estimator introduced in [S] . 

6.2 Indicator and Estimator Based on BDM 

For all F G £j and for i, j G {s, e}, let 

/^5,f = and I3ij,p = + fltj,p. 

and let 




it^ss,p 4~ flse,p) t^ee,F if^se,F l^ee,F^ Pse-a 

flss,pfdee,p -fl 

se,F 


and b, 


{(^se,p + flee,p) Pss,p — {(dss,p + flse,p) A 


se,F 


e,F 


flss,p/3ee,p fl. 


se,F 
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Using the basis functions ■0^ ^ and V’e f defined at the beginning of this section, a straightfor¬ 
ward calculation gives that 


/3: 


± 

SS,F 


^ - e,)"^ 


24|i^±| 


24\Kf\' 


and = 


- Sp) A (x^ -e^) 

48|K±| 


Solving the local problems in (4.13) with Vp = BDMp gives the following recovered flux in 
BDMg^Pf: 


^bdm _ 


E ^ a;\F\ct^-+ E (6.5) 


FeE, 


F&Er 




with the normal components of the recovered flux given by the weighted averages: 


= +(l-6.,^)d+ and = be,p + {1 - be,p) ap 


bdra 


The edge indicator has the following explicit formula: 

|d“-d+l |F|u;E F^£i, 

Mm ^ J 0, F e £d, 

\ 

,V2 


, I*; -S„,,I \F\ {p:.,, + •iP -,,+ PzA • 


where Wp is given by 


wp = {l- bA'^AA + 2(1 - &.,.)(! - be,c)/3“ + (1 - + bi p^t + .A 


Next, we present explicit formula of the element indicator in terms of the current 
approximation Up and geometrical information of elements. For each F C dK, denote by Fg 
and Fe the remaining two edges of K that is opposite to and e^, respectively. Then the 
indicator is computed by three terms as follows: 


^bdm ^ ^bdm\^^ ^ ^ VUp) ^ + (AVUp, VUp)^^^ . ( 6 . 6 ) 

The third term above is given in the previous section, and the other two terms may be computed 

by 

= Y. T. ir||r'l(B„,--D„, + M„,) 

FcdK F'CdK 

and {ap,Vup)p = - E TTlTU^r(xF)sign^(F) |F| (n^^L^,) . 

FcdK 
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Here, the and L^, have the following formulas: 


B 


FF> 


D 


FF' 


M 


FF' 


and 


p-bdm ^bdm ^ ^^F ’^f' ^ 

^S,F ^S,F' 48 1^1 

^bdm p^bdm ^^f ’^f' ^ 

^’F ^e,F, 48 1^1 

,^bdm phdnri ^ ^^F ’^f' 

^s,F^ ASiKl 


p^bdm ^bdm 

e,F e,Fi 48 |K| 


K'= E 

F'cdK 


t Fe > t ^ J sign^ (Fe) sign^ (Fg) I Fe 11 Fg I, 
sign;^(Fe) sign;^(Fg)|Fe| |Fg| 

- sign^(Fs)sign^(Fe)|Fs||Fe|, 

> ^Fi ) signx(Fg) IFs11F^|, 

signx(^e)l^e|t^, -<5-eJ;'signj^(Fg)|F,'|t^J . 


7 Numerical Experiments 


In this section, we report some numerical results for the Kellogg benchmark test problem HU. 
Let H = (—1,1)^ and 

u{r, 9) = r'^fj,{9) 


in the polar coordinates at the origin with fj,{9) being a smooth function of 9. The function 
u{r, 9) satisfies the diffusion equation in (2.1) with A = al, r^r = 0, / = 0, and 


a{x) 


R in (0,1)2 U (-1,0)2, 

1 in L!\([0,1]2 U [-1,0)2). 


The 7 depends on the size of the jump. In the test problem, 7 = 0.1 is chosen and is cor¬ 
responding to F « 161.4476387975881. Note that the solution u{r,9) is only in F^+'>'“'^(H) 
for some e > 0 and, hence, it is very singular for small 7 at the origin. This suggests that 
refinement is centered around the origin. 

This problem is tested by the standard ZZ estimator and its variation: 


= llVw. 7 - - G(Vri^)||o,n and = ||a^/ 2 v^^ _ q, ^/‘^G{aVu^)\\o^n- 


Here, is the standard ZZ estimator, i.e., the norm of the difference between the numerical 
and recovered gradients; and the is a modified version, where the flux is recovered in 
continuous finite element space. Both versions of the ZZ estimators perform badly with many 
unnecessary over-refinements along the interfaces (see Figures 3 and 4). 

Meshes generated by and for both RT and BDM recovery are shown in Figures 
iBa and , respectively. The refinements are centered at the origin and there is no over¬ 


refinements along the interfaces. Similar meshes for this test problem generated by other error 
estimators can be found in PEIE]. The comparisons between the tru e error in the energy 


norm and the estimators, ^ and are shown in Figures 10, and 12, respectively. All the 


estimators have effectivity indexes very close to one. Here, the effectivity index is defined as 
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Figure 3: mesh generated by the ZZ indicatot Figure 4: mesh generated by the modified ZZ 

indicator 



Energy error 



Figure 5: mesh generated by 


Figure 6: error and estimator 


the ratio of the estimator and the true error in the energy norm. Moreover, the slope of the 
log(dof)- log (the relative error) for both ^ and ^ are very close to —1/2, which indicates the 
optimal decay of the error with respect to the number of unknowns. 
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